library(sf) #for polygons
library(tidyverse) #for data wrangling
library(tmap) #for interactive maps
library(janitor) #for cleaning names
#create list of years
kern_field_years <- c(2016:2019)

#loop through read-in process
for (i in kern_field_years) {
  filename <- paste0("kernfields_", i)
  wd <- paste0("data/kernfields/kern", i, ".shp")
  assign(filename, read_sf(wd))
}

#select region (based off of farmer)
kern_subset <- kernfields_2019 %>% 
  filter(AGENT == "McMANUS/WILSON,MICHELE/AARIN")

#crop several other years
kern_2017_crop <- st_crop(kernfields_2017, kern_subset)

kern_2018_crop <- st_crop(kernfields_2018, kern_subset)

kern_2019_crop <- st_crop(kernfields_2019, kern_subset)

#create list of cropped fields and tidy
kern_cropped_list <- lapply(mget(sprintf("kern_%d_crop", 2017:2019)), 
                 function(x) 
                   {
                   clean_names(x) %>%
                     st_make_valid()
                     }
                 )

#bind cropped data together, include new column "source" saying what dataset it came from
sample_bind <- do.call(rbind, lapply(names(kern_cropped_list), 
                                    function(x) 
                                      cbind(kern_cropped_list[[x]], 
                                            source = x)
                                    )) %>% 
    #turn source column into year column by keeping only numbers
    mutate(year = gsub("\\D", "", source)) %>% 
    #Create ID column
    mutate(id = row_number()) %>% 
    #merge ID column to pmt_site
    mutate(pmt_site_id = str_c(pmt_site, "-", id)) %>% 
    #merge this column with year
    mutate(pmt_site_string = str_c(year, "_", pmt_site_id))



#select only column we need, clean geometry
sample_clean <- sample_bind %>% 
  dplyr::select(pmt_site_string) %>% 
  st_make_valid()
  

#write to shapefile
#st_write(sample_clean, "data/shapefiles_written/sample_clean.shp", append = F)

I took the shapefile sample_clean.shp and uploaded it to mapshaper’s online interface (https://mapshaper.org/). Then I entered in the console $ mosaic and downloaded results as a shapefile, as seen in the next code chunk.

This thread could potentially improve workflow https://github.com/mbloch/mapshaper/issues/353

#read in mapshaper output
sample_mosiac <- st_read("data/mapshaper_outputs/sample_clean/sample_clean.shp")
## Reading layer `sample_clean' from data source 
##   `/Users/irisfoxfoot/Desktop/Kern RA/kernRA_coding_projects/data/mapshaper_outputs/sample_clean/sample_clean.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 817 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 6159540 ymin: 2352351 xmax: 6187962 ymax: 2386672
## Projected CRS: NAD83 / California zone 5 (ftUS)
#intersect with field data 
sample_join <- st_intersection(st_buffer(sample_mosiac, dist = -1), sample_clean)

#pivot so each column is a year
sample_pivot <- sample_join %>% 
  separate(pmt_site_string, into = c("year", "pmt_site"), sep = "_") %>% 
  select(-FID) %>% 
  pivot_wider(names_from = year, values_from = pmt_site) %>%
  st_sf() %>% 
  st_buffer(dist = 1) %>%
  mutate(area_sq_ft = st_area(.))

#unnest listed pmt_site so they are each in own column
sample_pivot_sep <- sample_pivot %>% 
  unnest_wider(c(`2017`, `2018`, `2019`), names_sep = "_")

#write to csv
#write_csv(sample_pivot_sep, "data/shapefiles_written/2017_2019_sample.csv")

#write to shp file
#st_write(sample_pivot_sep, "data/shapefiles_written/2017_2019_shapefile.shp")

How to get rid of slivers?

Function to test how many true pmt_site_ids are dropped

pmt_deleted <- function (pivot_sep_name) {
  #create df of pivot table without geometry
  pivot_sep_no_geo <- pivot_sep_name %>% 
    select(-geometry, -area_sq_ft)
  
  #create function to check if pmt_site is present anywhere in the pivot table
  pmt_present <- function(x) 
    {
    present <- any(pivot_sep_no_geo==x)
    return(present)
  }
  
  #create list of original pmt_sites
  pmt_list <- as.list(sample_bind$pmt_site_id)

  #run through list and check if they are in pivot df
  present_list <- lapply(pmt_list, pmt_present)
  
  #create df of results
  present_df <- as.data.frame(cbind(pmt_list, present_list)) %>% 
    filter(is.na(present_list))

  return(present_df)
  
  }

Drop fields under 1 acre

#get rid of small fields
sample_pivot_over_acre <- sample_pivot %>% 
  filter(area_sq_ft >= units::set_units(43560,"ft^2"))

#unnest listed pmt_site so they are each in own column
sample_pivot_over_acre_sep <- sample_pivot_over_acre %>% 
  unnest_wider(c(`2017`, `2018`, `2019`), names_sep = "_")

#what fields are under one acre?
under_acre <- sample_bind %>%
  st_sf() %>% 
  mutate(area = st_area(.)) %>% 
  filter(area < units::set_units(43560,"ft^2"))

#gets rid of many but not all slivers
#deletes 65 true fields but, there are only 60 true fields that are under one acre. 
#what are the other five?

Negative buffer

#intersect with field data BUT negative buffer field data
neg_sample_join<- st_intersection(st_buffer(sample_mosiac, dist = -20), sample_clean)

#pivot so each column is a year
neg_sample_pivot <- neg_sample_join %>% 
  separate(pmt_site_string, into = c("year", "pmt_site"), sep = "_") %>% 
  select(-FID) %>% 
  pivot_wider(names_from = year, values_from = pmt_site) %>%
  st_sf() %>% 
  filter(!st_is_empty(.)) %>% 
  st_buffer(dist = 1) %>%
  mutate(area_sq_ft = st_area(.))

#unnest listed pmt_site so they are each in own column
neg_sample_pivot_sep <- neg_sample_pivot %>% 
  unnest_wider(c(`2017`, `2018`, `2019`), names_sep = "_")

#is every permit_site represented in the pivot table?

#with a buffer of -10 we miss 16 fields
#with a buffer of -20 we miss 31 fields
#with a buffer of -30 we miss 40 fields
#with a buffer of -50 we miss 76 fields

Conditional statement

#create column indicating if it was duplicated, area column
sample_conditional <- sample_join %>%
  group_by(pmt_site_string) %>% 
  mutate(dup = n()>1) %>% 
  ungroup() %>% 
  mutate(area_sq_ft = st_area(.))

#filter so fields smaller than one acre are only dropped if duplicated
sample_conditional_drop <- sample_conditional %>% 
  filter(!(dup  == "TRUE" & area_sq_ft <= units::set_units(43560,"ft^2")))

#pivot so each column is a year
conditional_sample_pivot <- sample_conditional_drop %>% 
  separate(pmt_site_string, into = c("year", "pmt_site"), sep = "_") %>% 
  select(-FID) %>% 
  pivot_wider(names_from = year, values_from = pmt_site) %>%
  st_sf() %>% 
  filter(!st_is_empty(.)) %>% 
  st_buffer(dist = 1) %>%
  mutate(area_sq_ft = st_area(.))

#unnest listed pmt_site so they are each in own column
conditional_sample_pivot_sep <- conditional_sample_pivot %>% 
  unnest_wider(c(`2017`, `2018`, `2019`), names_sep = "_")

#drops 25

Original Data

#reference what output should be
tmap_mode("view")

tm_shape(sample_bind) +
  tm_polygons(id = "pmt_site") +
  tm_facets(by = "year", sync = T) +
  tm_layout(panel.labels = c("2017", "2018", "2019"))

PMT_SITE Pivot Table

#test output
tmap_mode("view")

tm_shape(sample_pivot) +
  tm_polygons(id = "pmt_site_year") +
  tmap_options(check.and.fix = TRUE)

PMT_SITE Pivot Table (over 1 acre)

#test output
tmap_mode("view")

tm_shape(sample_pivot_over_acre) +
  tm_polygons(id = "pmt_site_year") +
  tmap_options(check.and.fix = TRUE)

PMT_SITE Pivot Table (conditionally over 1 acre)

#test output
tmap_mode("view")

tm_shape(conditional_sample_pivot) +
  tm_polygons(id = "pmt_site_year") +
  tmap_options(check.and.fix = TRUE)

PMT_SITE Pivot Table (negative 20 ft buffer)

#test output
tmap_mode("view")

tm_shape(neg_sample_pivot) +
  tm_polygons(id = "pmt_site_year") +
  tmap_options(check.and.fix = TRUE)